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Abstract 


Over  the  last  months,  much  effort  has  been  put  into  further  development  of  the  XBeach 
model.  A  visit  was  made  to  Delaware  by  Ap  van  Dongeren,  to  consult  with  Brad  Johnson  of 
ERDC  about  the  coupling  between  XBeach  and  Shorecirc;  it  was  decided  to  focus  on  further 
(joint)  development  of  XBeach.  The  following  improvements  and  extensions  were  carried 
out: 

•  Improvement  of  numerical  scheme  in  line  with  Stelling  and  Duinmeijer  method,  to 
improve  long-wave  runup  and  backwash  on  the  beach.  The  momentum-conserving 
form  is  applied,  while  retaining  the  simple  first-order  approach.  The  resulting 
scheme  has  been  tested  against  the  well-known  Carrier  and  Greenspan  test. 

•  Implementation  of  the  Generalised  Lagrangean  Mean  (GLM)  approach  to  represent 
the  depth-averaged  undertow  and  its  effect  on  bed  shear  stresses  and  sediment 
transport,  cf.  Reniers  et  al.  (2004) 

•  Implementation  of  Roelvink  (1993)  wave  dissipation  model  for  use  in  the 
instationary  wave  energy  balance  (in  other  words,  when  the  wave  energy  varies  on 
the  wave  group  timescale) 

•  Implementation  of  Soulsby  -  Van  Rijn  transport  formulations,  cf  Reniers  et  al 
(2004). 

•  Automatic  time  step  based  on  Courant  criterion,  with  output  at  fixed  time  intervals. 

•  Improvement  of  avalanching  mechanism,  with  separate  criteria  for  critical  slope  at 
wet  or  dry  points. 

A  first  series  of  (mainly  ID)  validation  tests  have  been  carried  out,  culminating  in  the 
successful  simulation  of  a  dune  erosion  test  carried  out  in  a  large-scale  wave  flume  (LIP11D 
test  2E,  Arcilla  et  al.,  1993).  The  results  are  briefly  shown  in  this  report. 


1  Introduction 


This  report  is  the  second  progress  report  of  the  project  ‘Modeling  of  Elurricane  Impacts’, 
contract  no.  N62558-06-C-2006,  which  was  granted  by  the  US  Army  Corps  of  Engineers, 
Engineer  Research  and  Development  Center  (ERDC),  European  Research  Office  and 
administered  by  FISC  SIGONELLA,  NAVAL  REGIONAL  CONTRACTING  DET 
LONDON,  SHORE/FLEET  TEAM. 

The  project  is  being  carried  out  by  Prof.  Dano  Roelvink  of  UNESCO-IHE  (Principal 
Investigator),  Dr.  Ad  Reniers  and  Jaap  van  Thiel  de  Vries  of  Delft  University  of  Technology 
and  Dr.  Ap  van  Dongeren,  Dirk-Jan  Walstra  and  Jamie  Lescinski  of  WL  |  Delft  Hydraulics. 

Over  the  last  months,  much  effort  has  been  put  into  further  development  of  the  XBeach 
model.  A  visit  was  made  to  Delaware  by  Ap  van  Dongeren,  to  consult  with  Brad  Johnson  of 
ERDC  about  the  coupling  between  XBeach  and  Shorecirc;  it  was  decided  to  focus  on  further 
(joint)  development  of  XBeach. 

A  number  of  improvements  and  extensions  were  made;  they  are  described  in  Chapter  3. 
Chapter  4  contains  an  update  of  all  formulations  used  and  their  numerical  schematization.  In 


1  -  1 


Netherlands  Centre  for  Coastal  Research 


Modeling  of  hurricane  impacts 
Progress  Report  1  -  March-May,  2006 


Chapter  5,  a  series  of  validation  tests  is  described,  from  analytical  experiments  to  check  the 
numerical  behaviour  to  a  large-scale  dune  erosion  experiment  in  a  wave  flume  (Arcilla  et  al, 
1993).  In  Chapter  6  we  draw  conclusions  and  sketch  our  plans  for  the  coming  period. 


2  Visit  Ap  van  Dongeren  to  Delaware 


In  a  two-day  visit  by  Ap  van  Dongeren  to  the  University  of  Delaware,  he  discussed  progress 
of  XBeach  and  ShoreCirc  with  Brad  Johnson.  Given  the  overlapping  functionality  in  the 
surf  zone  and  the  dedicated  functionality  of  XBeach  in  the  swash  zone  it  was  decided  to  not 
interface  with  Shorecirc,  but  rather  to  incoiporate  elements  of  Brad’s  code  into  XBeach.  To 
this  end,  the  source  code  of  XBeach  that  was  delivered  after  the  first  three  months  was 
discussed  in  detail,  including  more  recent  updates. 


3  XBeach  model  development  and 

validation 

3.1  Status  in  June  2006 

The  main  objective  of  the  XBeach  model  is  to  provide  a  robust  and  flexible  environment  in 
which  to  test  morphological  modelling  concepts  for  the  case  of  dune  erosion,  overwashing 
and  breaching.  The  top  priority  is  to  provide  numerical  stability;  first  order  accuracy  is 
accepted  since  there  is  a  need  for  small  space  steps  and  time  steps  anyway,  to  represent  the 
strong  gradients  in  space  and  time  in  the  nearshore  and  swash  zone.  Because  of  the  many 
shock-like  features  in  both  hydrodynamics  and  moiphodynamics  we  choose  upwind 
schematizations  as  a  means  to  avoid  numerical  oscillations  which  can  be  deadly  in  shallow 
areas. 

The  modelling  environment  should  be  flexible  and  the  code  easy  to  comprehend  and 
concise;  therefore  we  have  adapted  the  Matlab  environment  as  development  environment; 
conversion  to  Fortran  90/95  at  a  later  stage  will  be  straightforward. 

The  code  has  the  following  functionality: 

•  Depth- averaged  shallow  water  equations  including  time-varying  wave  forcing 
terms;  combination  of  sub-  and  supercritical  flows; 

•  Time- varying  wave  action  balance  including  refraction,  shoaling,  current  refraction 
and  wave  breaking; 

•  Wave  amplitude  effects  on  wave  celerity; 

•  Depth-averaged  advection-diffusion  equation  to  solve  suspended  transport; 

•  Bed  updating  algorithm  including  possibility  of  avalanching; 

•  Possibility  to  extend  to  parallel  multi-domain  version; 
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3.2  Improvements  and  extensions  June-August  2006 

The  following  improvements  and  additions  have  been  implemented  and  tested  in  the  last 
three  months: 

•  Improvement  of  numerical  scheme  in  line  with  Stelling  and  Duinmeijer  method,  to 
improve  long-wave  runup  and  backwash  on  the  beach.  The  momentum-conserving 
form  is  applied,  while  retaining  the  simple  first-order  approach.  The  resulting 
scheme  has  been  tested  against  the  well-known  Carrier  and  Greenspan  test. 

•  Implementation  of  the  Generalised  Lagrangean  Mean  (GLM)  approach  to  represent 
the  depth-averaged  undertow  and  its  effect  on  bed  shear  stresses  and  sediment 
transport,  cf.  Reniers  et  al.  (2004) 

•  Implementation  of  Roelvink  (1993)  wave  dissipation  model  for  use  in  the 
instationary  wave  energy  balance  (in  other  words,  when  the  wave  energy  varies  on 
the  wave  group  timescale) 

•  Implementation  of  Soulsby  -  Van  Rijn  transport  formulations,  cf  Reniers  et  al 
(2004). 

•  Automatic  time  step  based  on  Courant  criterion,  with  output  at  fixed  time  intervals. 

•  Improvement  of  avalanching  mechanism,  with  separate  criteria  for  critical  slope  at 
wet  or  dry  points. 

3.3  General  setup  (update) 

The  program  XBeach  consists  of  a  main  Matlab  script,  xbeach.m,  and  a  number  of functions 
that  operate  on  two  structures : 

•  par  -  this  contains  general  input  parameters 

•  s  -  this  contains  all  the  arrays  for  a  given  computational  domain 

For  a  single-domain  run,  one  structure  s  is  passed  between  flow,  wave,  sediment  and  bed 
update  solvers,  which  extract  the  arrays  they  need  from  the  structure  elements  to  local 
variables,  do  their  thing  and  pass  the  results  back  to  the  relevant  structure  elements.  This 
makes  the  overall  program  clear,  prevents  long  parameter  lists  and  makes  it  easy  to  add 
input  variables  or  arrays  where  needed. 

For  multi-domain  runs,  one  can  define  multiple  instantiations  of  the  structure  s  which  are 
passed  to  the  same  functions;  an  additional  function  is  needed  to  pass  the  boundary 
information  between  the  domains  back  and  forth.  We  have  carried  out  a  simple  test  of  this 
principle,  without  actually  implementing  a  multi-processor  version,  which  confirms  that  the 
data  structure  can  handle  this  case. 

In  the  Table  1  we  will  outline  the  various  functions  and  their  purposes.  The  main  script 
xbeach.m  is  reproduced  in  Table  2  below. 


Function  call 

Purpose 

[par]  =  waveinput 

Creates  elements  of  structure  par  containing  wave  input 
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parameters 

[par]  =  flowinput(par) 

Adds  elements  of  structure  par  containing  flow  input 
parameters 

[par]  =  sedinput(par) 

Adds  elements  of  structure  par  containing  sediment  input 
parameters 

[s]  =  gridbathy; 

Creates  grid  and  bathymetry  and  stores  them  in  structure  s 

[s]  =  wave_dist  (s,par); 

Creates  initial  directional  spectrum  at  sea  boundary 

[s,par]=wave_init  (s,par); 

Initialises  arrays  (elements  of  s)  for  wave  computations 

[s]  =  flow  init  (s,par); 

Initialises  arrays  (elements  of  s)  for  flow  computations 

[s]  =  sed_init  (s); 

Initialises  arrays  (elements  of  s)  for  sediment  computations 

[s]  =  wave_bc  (s,par,it); 

Wave  boundary  conditions  update,  each  timestep 

[s]  =  flow_bc  (s,par,it); 

Flow  boundary  conditions  update,  each  timestep 

[s]=wave_timestep  (s,par); 

Carries  out  one  wave  timestep 

[s]=flow_timestep  (s,par); 

Carries  out  one  flow  timestep 

[s]=transus(s,par); 

Carries  out  one  suspended  transport  timestep 

[s]=bed_update(s,par); 

Carries  out  one  bed  level  update  timestep 

output(it,s,par) 

Performs  output 

Table  1  Overview  of  Matlab  functions  XBeach 


Of  these  functions,  the  input  functions  and  grid  bathy  are  case-dependent  at  present  and 
actually  contain  the  input  values  and  grid  and  bathymetry  definitions.  This  can  however 
easily  be  replaced  by  functions  with  the  same  output,  which  read  data  from  input  files  or 
input  them  through  a  screen  dialog.  We  aim  to  leave  the  initialisation,  boundary  conditions 
and  computational  functions  case-independent.  The  output  function  can  be  adapted  to  fit 
specific  needs. 
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clear  all 

"o 

%  General  input  per  module 
[par]  =  wave_input; 

[par]  =  f low_input (par ) ; 

[par]  =  sed_input (par )  ; 

%  Grid  and  bathymetry 
[s]  =  grid_bathy; 

%  Directional  distribution  wave  energy 
[s]  =  wave_dist  (s,par); 

%  Initialisations 

[s,par]  =  wave_init  (s,par); 

[s]  =  flow_init  (s,par); 

[s,par]  =  sed_init  (s,par); 

nt=par . nt ; 

par . tnext=par . tint; 

it=0; 

while  par . t<par . tstop; 

%  Wave  boundary  conditions 
[par,s]  =  wave^bc  (s, par, it); 

%  Flow  boundary  conditions 
[s]  =  flow_bc  (s, par, it); 

%  Wave  timestep 

[s]  =  wave_timestep  (s,par); 

%  Flow  timestep 

[s,par]  =  f low^timestep  (s,par); 

%  Suspended  transport 
[ s ] =transus ( s , par ) ; 

%  Bed  level  update 
[s] =bed_update (s,par)  ; 

%  Output 

[ it , s ] =output ( it, s, par ) ; 

End 

Table  2  Main  routine  xbeach.m 


3.4  Overview  of  modifications  Matlab  code 


Routine 

Modifications  June-August  2006  (from  vOOO  to  vOlO) 

Baldock 

None 

Bedupdate 

Avalanching 

Morphological  factor 

Dispersion 

None 

Flowbc 

Timestep  management 

Flow  int 

Include  initial  water  level 

Initialize  qx,qy 

Flowinput 

Added  initial  water  level,  Tstart,  Tint,  Tstop 

Flowtimestep 

GLM  velocities 

Stelling  and  Duinmeijer  scheme  (momentum-conserving,  first  order) 
Bug-fix  seaward  boundary 

Automatic  timestep 

Grid  bathy 

Model  dependent 

Output 

Model  dependent 

Roelvink 

New,  computes  dissipation  according  to  Roelvink  (1993) 
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Sbvr 

New;  port  formulations  according  to  Soulsby-  van  Rijn 

Sed_init 

None 

Sedinput 

4  parameters  added 

Slope 

None 

Transus 

Eulerian  velocities  used 

Equil.  Concentration  according  to  Soulsby-  van  Rijn 

Wavebc 

Model-dependent  treatment  of  wave  boundary  conditions  (e.g.  energy 
time-series  read  from  file) 

Wave  init 

Intriduced  initial  water  level 

Set  initial  value  time  step 

Wave  input 

Added  parameters  Roelvink  ’93  formulation 

Added  choice  of  formulation 

Wave_timestep 

Added  Roelvink  ’93  formulation 

Computation  of  Urms,  Ustokes 

Xbeach 

Time  management  introduced 

3.5  Conversion  to  Fortran  90/95 

The  reason  to  work  on  a  Fortran  version  of  the  code  is  that  it  is  easier  to  create  stand-alone 
versions  on  any  platform,  the  code  may  be  easier  to  parallellize  and  it  may  be  faster. 

Jamie  Lescinski  at  Delft  Flydraulics  has  converted  the  Matlab  code  to  Fortran90.  She  is 
trying  two  tracks: 

•  to  convert  the  Fortran  code  to  MEX  files  and  call  them  from  a  Matlab  environment, 
This  is  posing  some  problems  that  may  be  due  to  a  bug  in  Matlab,  which  she  is 
presently  trying  to  solve  with  the  Matlab  help  desk; 

•  to  convert  all  code  to  Fortran.  This  has  just  been  completed  and  is  ready  to  be 
tested. 

The  conversion  from  Matlab  to  Fortran  appears  to  be  relatively  easy,  as  the  ‘derived  type’  in 
Fortran  90/95  is  very  similar  to  Matlab  structures. 
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4  XBeach  formulations  (update) 


4.1  Wave  action  equation  solver 

The  wave  forcing  in  the  shallow  water  momentum  equation  is  obtained  from  a  time 
dependent  version  of  the  wave  action  balance  equation.  Similar  to  Delft  University’s 
HISWA  model,  the  directional  distribution  of  the  action  density  is  taken  into  account 
whereas  the  frequency  spectrum  is  represented  by  a  single  mean  frequency.  The  wave  action 
balance  is  then  given  by: 


dA  dcA  dc  A  dcnA  D 

- 1 - - — i - - — i - - —  = - 

dt  dx  dy  dO  a 


with  the  wave  action: 


A(x,y,0 ) 


Sw(x,y,0) 

<r(x,y) 


(4.1) 


(4.2) 


where  Sw  represents  the  wave  energy  in  each  directional  bin  and  a  the  intrinsic  wave 
frequency.  The  wave  action  propagation  speeds  in  x-  and  y-direction  are  given  by: 


cx(x,  y,  0)  =  cg(x,  v> cos(#)  +  u(x,  y) 
Cy  (x,y,0)  =  cg  (x,  yy  sin(#)  +  v(x,  y) 


where  6  represents  the  angle  of  incidence  with  respect  to  the  x-axis.  The  propagation  speed 
in  6  -space  is  obtained  from: 


cg(x,y,0)  = 


a 


sinh  2kh 


8h  . 


dh 


—  sin  6 - cos  6 

dx  dy  j 


f 


+  cos  6 


,du 


du 


sin  6 - cos  6  — 


f 


TSUI# 


dx 
,  dv 


dy 


+ 


•  ndV 

sin# - cos#  — 


(4.4) 


dx 


dy 


taking  into  account  bottom  refraction  (first  term  on  the  RHS)  and  current  refraction  (last  two 
terms  on  the  RHS).  The  wave  number  k  is  obtained  from  the  eikonal  equations: 


dk  dco 
— ^  +  —  =  0 
dt  dx 

dky  !  da)  _  Q 

dt  dy 


(4.5) 
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where  the  subscripts  refer  to  the  direction  of  the  wave  vector  components  and  CD  represents 
the  absolute  radial  frequency.  The  wave  number  is  the  obtained  from: 


k  ~  yjkx  +  k 

The  absolute  radial  frequency  is  given  by: 


co  =  <j  +  k»u 

and  the  intrinsic  frequency  is  obtained  from  the  linear  dispersion  relation: 


(4.6) 


u  =  yjgk  tanh  kh 


(4.7) 

(4.8) 


The  group  velocity  is  obtained  from  linear  wave  theory: 

r\  kh  ' 


cg  =  nc  = 


- b 

y2  sinh2  kh 


a 

k 


(4.9) 


This  concludes  the  advection  of  wave  action.  The  wave  energy  dissipation  due  to  wave 
breaking  is  modelled  according  to  Baldock  et  al.  [1998]: 


D  =  -aQtpgf,{Hl  +  H L) 


(4.10) 


with  a  =  0(1)  and  fm  representing  the  mean  intrinsic  frequency.  The  fraction  of  breaking 
waves  is  given  by: 


Qh  =  exp 


f  hi  y 

TT 2 

\  1  rms  J 

(4.11) 


where  the  breaking  wave  height  is: 


rr  0-88  , 

H,  = - tanh 


ykh 

088 


(4.12) 


and  y  is  a  calibration  parameter.  The  root  mean  square  wave  height  is  obtained  from: 


H 

1  rms 


8fSw(x,y,0)d0 


Pg 


Pg 


(4.13) 


Next  the  total  wave  dissipation,  D  ,  is  distributed  proportionally  over  the  wave  directions: 

(4.14) 


EJx,y) 
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This  closes  the  set  of  equations  for  the  wave  action  balance.  Given  the  spatial  distribution  of 
the  wave  action  and  therefore  wave  energy  the  wave  forcing  can  be  calculated  utilizing  the 
radiation  stress  tensor: 


And: 


r  dS 

dS  ] 

F=- 

XX  | 

xy 

X 

v  dx 

dy  J 

(  dSxv 

esm ) 

F=-\ 

y  1 

yy 

y  \ 

v  dx 

dy  a 

ff—  (l  +  cos^  6 

>)-±l 

/  0 

1 ) 

^=^  =  Jsin6>c°s6» 


sje 

\ 

c 


vie 


S;),=j(A(1  +  siirP)-b 


SjdO 


(4.15) 


(4.16) 


Similar  to  the  solution  of  the  shallow  water  equations  we  use  an  up-wind  schematisation  to 
solve  the  wave  action  balance.  The  wave  action  is  given  at  the  same  points  at  the  water 
level.  The  advection  of  wave  action  is  then  discretized  as  follows: 


An  rn  A”  _  rn  A" 

OCxA  r,  ;  w  _  Lx,i,j,k^iJ,k 


dx 


(i,j,k) 


,<^>0 


a  »  A"  r"  A"  -  rn  A" 

OCx A  I'j  j  |.  j  _  x,i+\,j,kn-i+\,j,k  cxj,j,k‘ni,j,k  C>1  <  q 

dx  ”  Xmj-Xij 


x,i,j,k 


(4.17) 


dcnA" 

y 

dy 


(i,j,k)  = 


rn  An  _ rn  An 

y,i,j,k  i,j,k  y  ,k  ij -1  ,k 


A"  -rn  A” 

k  yl;,/+U  ^y,i,i,k  i,i 


o 


r)rn  A "  r"  - 

y  - (i,j ,k)  =  y^k^-j+ !■*  ^yMX^iJJk  ? cn  (  <() 


dy 


Tij+i  -  Tij 


'  y,i,j,k 


(4.18) 


a"  rn  A"  —  r "  A" 

OCgJi  /.  .  _  ^e,i,j,k  i,Uk  L'dJJ,k-l  i,j,k-\  „n 


do 


(i,j,k): 


^ij.k 


’  C9,i,j,k  >  ^ 


a  »  A"  r"  A"  —  r"  A " 

=  0JJM1  ,-,M+1 - e,,J,k  ,J,k  n  0 


(4.19) 


Similar  for  the  wave  action  balance: 


A"+ 1  _  An 

At 


dcnxA"  dc;A"  dc"eAn  _D 
dx  ij,k  dy  i  jk  dO  tj,k  <Jijx 


(4.20) 
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which  yields  the  wave  energy  at  the  new  time  level. 


4.2  Shallow  water  equations  solver 


Shallow  water  equations,  neglecting  Coriolis  and  horizontal  diffusion  terms,  and  (grey 
terms),  for  the  moment,  wind  shear  stress: 


du  du  du  t.  drj  Fr 

—  +  u  —  +  v — =  — - —~g — -  +  — 

dt  5x  dy  ph  ph  dx  ph 

(4.21) 

dv  dv  dv  Tsy  Tbv  drj  Fy 

- h  U - h  V - —  H - - g - 1 - 

dt  dx  dy  ph  ph  dy  ph 

(4.22) 

dn  dhu  dhv 

—  + - + - =  0 

dt  dx  dy 

(4.23) 

Here,  h  is  the  water  depth,  u,  v  are  velocities  in  x  and  and y  direction,  Tbx,Tbv  are  the  bed 
shear  stresses,  g  is  the  acceleration  of  gravity,  lj  is  the  water  level  and  Fx,Fv  are  the 
wave-induced  stresses. 


We  apply  an  upwind  schematisation,  since  the  horizontal  scale  of  the  problem  is  limited  and 
such  a  scheme  deals  with  shocks  in  a  natural  way. 

We  apply  a  staggered  grid,  where  bed  levels  and  water  levels  are  defined  in  the  centre  of 
cells,  and  velocity  components  at  the  cell  interfaces. 

If  nx,ny  are  the  number  of  cells  in  both  directions,  the  water  level  points  are  numbered  from 
1  to  nx+1  and  from  1  to  ny+1. 

The  water  level  gradients  are  computed  at  the  cell  interfaces  and  are  given  by: 


fh  y)  =  Mzi 

OX 


dy 


(4.24) 


(4.25) 


For  computing  the  shear  stresses  at  the  cell  interfaces  we  need  the  velocity  magnitudes  at 
these  interfaces.  These  are  composed  by  combining  the  normal  velocity  component  at  the 
interface  and  the  average  of  the  4  adjacent  tangential  components: 

1 

lv  +vij  ■  •/.!./  i  ■  ri+lJJ 

(4.26) 


vu,i,j  =  +  vi,j  +  VMJ-1  +  vm  j) 


uv,ij  =  +  uu  +  i  +  UiJ+ 1) 


The  water  depth  in  each  cell  is  computed  as: 

Kj=r!ij-Zb,i,j 


(4.27) 
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For  the  depth  at  cell  interfaces,  following  Stelling  and  Duinmeijer  (2003)  we  distinguish 
between  the  depth  used  in  the  continuity  equation  and  that  used  in  the  momentum  equation. 
The  depth  at  the  interfaces  for  the  continuity  equation  is  taken  as  the  upwind  depth  in  case 
the  velocity  is  greater  than  a  minimum  velocity,  or  the  average  depth  between  the  adjacent 
cells  in  case  the  velocity  is  less  than  this  minimum  velocity: 


Kj 

’  M/,y 

>  u 

min 

Kuj 

’  Uu 

<  —u 

mm 

AU.  \  <U  ■ 

9|  i,j\  min 

Kj 

vu 

>Vmin 

hj+ 1 

<  -Vmin 

KiJ=\(Kj+Kj+ 1)>  |Vh/|<Vmi„ 


(4.28) 


(4.29) 


For  the  depth  in  the  momentum  balance  we  take  the  average  depth  between  the  cell  centers: 


Ku,iJ=^(hiJ+hM,j) 

hmVjJ=^(hiJ+hj+ 1)’ 


(4.30) 

(4.31) 


The  advection  terms  in  x-direction  are  approximated  as  follows: 


du'  _  1  K.i.f,,  •  hu.i  I,/  Kj-f-K 

dxij  2 


mu,i,j 


n  n 


du"  _  1  Kj.,11,.,  •  »,  i./  <i j~Kj 

dxu  2  h . .  , 


xmj~xu 


Xj>  0 


Kj<  o 


(4.32) 


^  n  n  n 

du_  =  n  uiJ+1  - 

^  u,i,j  n  n 

i,j  y i,j+ 1  y i,j~ i 


5v 


The  advection  terms  in  y-direction  are  approximated  as  follows: 


(4.33) 


<9v  _  i  Kjjvij  +  vij  ~  vij- 1 


-,in  n  .  j  n  n  n  n 

dv  _  1  hv4JvUj  +  hVJJ+lviJ+1  v,  J+1  -  v;.  y 


hn  . . 

mv,i,j 


y,- 


j+ 1  -y 


’  vu  >  o 


,<.<0 


(4.34) 


3v  „  v,.+1  -  vi+1 

W -  =  U  .  . - - - - 

dx  V’,J  ~n 


-  ‘-j 


XMJ  T.l./ 


(4.35) 


The  momentum  equation  is  discretized  as  follows: 
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n+ 1  n 
Ui.j  -Ujj 


At 


/*,  n  ^  n 

du  du 

-u —  -v  — 

dxij  dry  Uj 


§u, 


n  n 


2  .  n  2 
J  +Vu,i,J 


KjjC 


n"  -  rf  F 
g-^lA — Fl  +  ^AL_  (4.36) 

pk,ij 


n+ 1  n 

V  —V 

FJ  l,J 


At 


n  ^  n 

OV  OV 

=  -v  —  -  u  — 
dytj  dxij 


gv. 


n  /  n 


2  ,  n  2 

+  vu 


h"  C 

nv,i,F 


g  //li1  'h  ■  FvJJ  (4.37) 
yoH-yij  Phij 


From  this,  the  velocities  at  the  new  time  step  level  are  computed.  The  water  level  is  then 
updated  by: 


r,''+1_r,«  Un+lhn  —  ;/"+1  h"  v',+1h  -v"+1/7" 

'li.j  '/i,j  _  UiJ  ,li,j  Ui-l,jni-\,J  Vl,j  nij  Vi,j-\ni,j-\ 


At 


•Vu,.i  •V,.i-i.i 


y v,ij  y v,i,j-i 


(4.38) 


Generalized  Lagrangian  Mean  formulation 

To  account  for  the  wave  induced  mass-flux  and  the  subsequent  (return)  flow  the  shallow 
water  equations  are  cast  into  a  Generalized  Lagrangian  Mean  (GLM)  formulation  (Walstra 
et  al,  2000).  To  that  end  the  Eulerian  shallow  water  velocity  uE  is  replaced  with  its 
lagrangian  equivalent,  u  : 


L  E  .  S  7  L  E  ,  S 

u  =u  +u  ana  v  =v  +v 


(4.39) 


and  us ,  v5  represents  the  Stokes  drift  in  x-  and  y-direction  respectively  (Phillips,  1977): 


uS  =  E^cose  and  vS  =  £> nP 

phc  phc 


(4.40) 


where  the  wave-group  varying  short  wave  energy  and  direction  are  obtained  from  the  wave- 
action  balance.  The  resulting  GLM-momentum  equations  are  given  by: 


duL  L  duL  i  duL 

- +  u  - +  V  - 

dt  dx  dy 

dvL  L  dvL  L  dvL 

- +  u  - +  V  - 

dt  dx  dy 


ph  dx  ph 
ph  dy  ph 


(4.41) 


for  the  x-  and  y-direction  respectively.  This  operation  shows  that  the  GLM  equations  for  the 
depth- averaged  flow  are  very  similar  to  the  previously  described  Eulerian  formulation,  with 
the  exception  of  the  bottom  shear  stress  terms  that  are  calculated  with  the  Eulerian  velocities 
as  experienced  by  the  bed: 


ELS  j  E  L  S 

u  =u  -u  and  v  =v  -v 


(4.42) 
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and  not  with  the  GLM  velocities.  Also,  the  boundary  condition  for  the  flow  computations 
has  to  be  expressed  in  functions  of  (uL  ,  vL  )  and  not  (uE,  vE). 


4.3  Sediment  transport 


Advection-diffusion  scheme 

The  sediment  transport  is  modeled  with  a  depth-averaged  advection  diffusion  equation 
[Gallapatti,  1983]: 


dhC  dhCuE 

- 1 - 

dt  dx 


dhCvE  d 

+ - +  — 

dy  dx 


Dhh 

h  dx 


(4.43) 


where  C  represents  the  depth-averaged  sediment  concentration  which  varies  on  the 
inffagravity  time  scale.  The  entrainment  of  the  sediment  is  represented  by  an  adaptation 
time  T„  given  by  a  simple  approximation  based  on  the  local  water  depth,  h,  and  sediment 
fall  velocity  ws: 


( 


T  =  max 


0.05  —  ,0.2 

w. 


J 


(4.44) 


where  a  small  value  of  T  corresponds  to  nearly  instantaneous  sediment  response.  The 
entrainment  or  deposition  of  sediment  is  determined  by  the  mismatch  between  the  actual 
sediment  concentration  and  the  equilibrium  concentration,  Ceq,  thus  representing  the  source 
term  in  the  sediment  transport  equation. 


The  differential  equations  for  the  advection  diffusion  of  sediment  is  solved  with  finite 
differences  using  the  first  order  up-wind  scheme  discussed  earlier  with  the  water  depths  at 
the  old  time  level  and  the  corresponding  velocities  at  the  new  time  level.  The  horizontal  x- 
advection  is  then  given  by: 


dhCuE  ^ 

i 

(  j  n /'-»«  E,n+ 1  ^ 

\htu  J 

i  (in  /-in  E,n+ 1  \ 

|  -\hCu 

i,j  v  'i-l  ,j 

l  dx  J 

Uj 

X 

:i.j~Xi  G 

f  dhCuE ) 

,  i 

[hnCnuE’n+l] 

| .  .-(A"CV"+1) 

V  dx  J 


IJ 


..Eji+1  .  r\ 

,  ui  j  >  0 


, uf  ’j  <  0  (4.45) 


a  similar  expression  for  the  horizontal  advection  in  the  y-direction: 
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fdhCvE A 


dy 


fdhCvE A 


<9y 


( h"C"vE’"+l ) 

|  -(h"CnvE’"+1) 

i,j  V  //-] 

y  i,j — y  ij-\ 

( hnCnvE’n+l ) 

)jj+i-(hnC"vE’n+1) 

hj 


'hj 


yij+i-yt. 


,vf;"+1  >0 


,  <  0  (4.46) 


The  horizontal  diffusion  is  evaluated  at  the  old  time  level  n  and  approximated  by: 


SCYl  {DHhCdx)  ~(DHhCdx). 


\dXy 


DHh 


dx 


JJ 


IJ 


XMJ  ~Xuj 


Where  the  cross-shore  gradient  in  the  sediment  concentration  is  given  by: 

f  dC^ 


C  -C 

_  J  uj 
\dx)iJ  XM  J~XiJ 


(4.47) 


(4.48) 


And  similarly  for  the  y-direction: 


r  a  c  sc 


8_ 


DHh 


dy 


W 


JJt 


(DHhCdy)ij+\  {DHhCdy)i] 

y^-yu 


’vo<0  (4.49) 


Where  the  along-shore  gradient  in  the  sediment  concentration,  Cy  ,  is  given  by: 


r  dC ^ 


C  -C 

_  1  ^ij 


\  dy  Jij  ytj+ 1  -  ytJ 

The  time  up-date  of  the  sediment  concentration  is  then  given  by: 

dhCuE  1  f  dhCvE 


hn+xCn+l-hn  C". 

uj  uj  UJ  UJ 


+ 


+ 


+ 


+ 


At 

L 

dx 

T 

L  dy 

1., 

~  d 

\DhKi 

-  n 

+ 

-  Uj 

"  d 

\Dh^  T 

■  n 

~hCeq-hC 

dx 

l  dx\ 

dy 

L  5yJ_ 

■  Uj 

L  u  J 

(4.50) 


(4.51) 


The  bed-update  is  discussed  next.  Based  on  the  gradients  in  the  sediment  transport  the  bed 
level  changes  according  to: 


dzh  dS ..  dSv 


dt  dx  dy 


=  0 


(4.52) 


where  Sx  and  Sv  represent  the  sediment  transport  rates  in  x-  and  y-direction  respectively, 
given  by: 
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and 


S*xjj  = 


rt/7  _ 

‘Vj  - 


dhCuL 

dx 


8hCv£ 


+ 


-‘i.j  L- 


dy 


+ 


-‘hj 


d_ 

dx 

8_ 

dy 


Dhh  — 
*  dx 


Dhh  — 
dy 


-»,J 


(4.53) 


(4.54) 


The  bed-update  is  then  approximated  by: 


n+ 1  n 
Zb,i,j  ~  Zb,i,j 


At 


+  f 

J  m 


rin  rin  rin 

^x,ij  -  ^ xj-lj  .  ^y,i,j  ~  ^yjj- 1 


Ax 


+ 


Ay 


=  0 


(4.55) 


where  fmor  represents  a  morphological  factor  to  speed  up  the  bed  evolution. 


Transport  formulations 

The  equilibrium  sediment  concentration  can  be  calculated  with  various  sediment  transport 
formulae.  At  the  moment  the  sediment  transport  formulation  of  Soulsby-van  Rijn  (Soulsby, 
1997)  has  been  implemented.  The  Ceq  is  then  given  by  : 


r  _  At  + 

h 


E  |2 


+0.018- 


c 


d  J 


o.5  y 

—  u 

cr 

J 


2.4 


(1  -  a,m) 


(4.56) 


where  sediment  is  stirred  by  the  Eulerian  mean  and  infragravity  velocity  in  combination 
with  the  near  bed  short  wave  orbital  velocity  obtained  from  the  wave-group  varying  wave 
energy.  The  combined  mean/infrgarvity  and  orbital  velocity  have  to  exceed  a  threshold 
value,  ucr,  before  sediment  is  set  in  motion.  The  drag  coefficient,  C*  is  due  to  flow  velocity 
only  (ignoring  short  wave  effects).  To  account  for  bed-slope  effects  on  the  equilibrium 
sediment  concentration  a  bed-slope  correction  factor  is  introduced,  where  the  bed-slope  is 
denoted  by  m  and  ah  represents  a  calibration  factor.  The  bed  load  coefficients  Asb  and  the 

suspended  load  coefficient  Ass  are  functions  of  the  sediment  grain  size,  relative  density  of 
the  sediment  and  the  local  water  depth  (see  Soulsby  [1997]  for  details). 


4.4  Bottom  updating 


Avalanching 

To  account  for  the  slumping  of  sandy  material  during  storm-induced  dune  erosion 
avalanching  is  introduced  to  update  the  bed-evolution.  Avalanching  is  introduced  when  a 
critical  bed-slope  is  exceeded: 
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dz. 


dx 


>  m. 


(4.57) 


Where  the  estimated  bed  slope  is  given  by: 


d ]Zb  _  Zb,M,j  Zb,i,j 


dx  Ax 

The  bed-change  within  one  time  step  is  then  given  by: 


(4.58) 
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(4.59) 


where  a  threshold  of  0.005  m  has  been  introduced  to  prevent  the  generation  of  large 
shockwaves.  The  corresponding  bed  update  is  given  by: 


7n+ 1 
' b.i,j 
n+ 1 


+  Az 


b,i,j 


=  Zn  —  At 

b,i+l,j  b,i+l,j  b,i,j 


(4.60) 


To  account  for  continuity,  e.g.  when  sand  is  deposited  within  the  wet  part  of  the  domain,  the 
water  level  is  also  updated: 


_rn+ 1 
n+ 1 


=  Z„ 


+  Az 


b,i,j 


z  ‘  ‘  =  z'1  —  Az 

^s,i+\J  s,i+\,j  b,i,j 


(4.61) 


Similar  expressions  are  used  for  the  subsequent  avalanching  in  the  y-direction. 
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5  XBeach  validation 


5.1  Long  wave  propagation  and  numerical  damping 

The  purpose  of  the  first  test  is  to  check  if  the  scheme  is  not  too  dissipative  and  that  it  does 
not  create  large  errors  in  propagation  speed. 

A  long  wave  with  a  small  amplitude  of  0.01  m  and  period  of  80  s  was  sent  into  a  domain  of 
5  m  depth,  grid  size  of  5  m  and  a  length  of  1  km.  At  the  end,  a  fully  reflecting  wall  is 
imposed.  The  wave  length  in  this  case  should  be  sqrt(9.81*5)*80  =  560  m.  The  velocity 
amplitude  should  be  sqrt(g/h)*amp  =  sqrt(9.81/5)*0.01  =  0.014  m.  After  the  wave  has 
reached  the  wall,  a  standing  wave  with  double  amplitude  should  be  created. 

As  Figures  5.1  and  5.2  show,  the  model  accurately  represents  this  situation.  There  is  hardly 
any  dissipation,  the  wave  length  is  very  close  to  what  it  should  be  and  there  is  no  re- 
reflection  off  the  seaward  boundary. 
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Figure  5.1  Snapshots  of  water  level  and  velocity  at  T/4  intervals;  just  as  wave  hits  wall. 


Netherlands  Centre  for  Coastal  Research 


5-17 


Modeling  of  hurricane  impacts 
Progress  Report  1  -  March-May,  2006 


Figure  5.2  Same  as  5.1,  after  long  time. 
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5.2  Long  wave  runup  on  sloping  beach;  comparison 
with  Carrier  and  Greenspan  (1  958) 

The  purpose  of  this  test  is  to  check  the  ability  of  the  model  to  represent  runup  and  rundown 
of  long  waves.  To  this  end,  a  comparison  was  made  with  the  analytical  solution  by  Carrier 
and  Greenspan  (1958),  which  describes  the  motion  of  harmonic,  non-breaking  long  waves 
on  a  plane  sloping  beach  without  friction. 

In  Figure  5.3  the  model  results  for  waves  at  an  amplitude  of  Vi  the  breaking  wave  amplitude 
are  shown,  at  1/20  T  intervals.  The  agreement  is  quite  good,  though  there  is  very  small 
disturbance/lag  during  rundown.  Typical  of  the  solution  is  that  the  profiles  during  runup  and 
rundown  should  be  exactly  equal;  apart  from  a  small  area  near  the  water  line  during 
rundown,  this  is  the  case. 


surface  elevation 


cross-shore  distance  (m) 


Figure  5.3  Surface  elevation  snapshots  at  1/20  T  intervals;  model  (thin  brown  lines)  and  analytical  solution  (thick 
blue  lines). 
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5.3  Stationary  wave  propagation,  dissipation  and 
setup 

The  purpose  of  this  test  was  to  check  the  wave  energy  and  momentum  balance  in  stationary 
mode. 

The  Delta  Flume  test  of  Arcilla  et  al.  (1993),  test  2E  was  used.  This  test  was  carried  out  with 
an  increased  water  level  and  significant  dune  erosion  occurred  during  the  test.  Besides, 
extensive  hydrodynamic,  sediment  transport  and  morphological  measurements  were  carried 
out. 

The  dissipation  model  of  Baldock  (1998)  was  applied,  with  a  gamma  value  of  0.8.  The 
results  show  a  good  agreement  for  the  Hrms  wave  height,  mean  setup  and  rms  value  of  the 
orbital  velocity. 
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Figure  5.4  Wave  height  decay,  setup  and  orbital  velocity,  modelled  in  stationary  mode  (drawn  blue  line)  vs. 
measured  (red  dots).  LIP  1  ID  test  2B. 
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5.4  Instationary  surf  zone  flows  in  large-scale  flume 
test 

The  purpose  of  this  test  was  to  verify  the  hydrodynamics  of  the  model  when  run  in 
instationary  mode,  viz.  with  a  time-varying  wave  energy  imposed  at  the  offshore  boundary. 

The  same  test  case  as  before  was  used.  The  time  series  of  wave  energy  was  generated  by  a 
procedure  described  in  Roelvink  (1993b),  based  on  a  JONSWAP  spectral  shape  as  was 
applied  in  the  test.  Zero-order  steering  was  applied  in  the  flume  test;  therefore  no  incident 
bound  long  wave  was  imposed  in  the  numerical  experiment. 

The  results  are  given  in  Figures  5.5  and  5.6,  and  show  that  the  short  wave  decay  and  the 
generated  long  waves  are  represented  quite  well,  both  in  surface  elevation  and  in  near¬ 
bottom  velocity.  Also  the  time-averaged  water  level  is  represented  quite  accurately. 
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Figure  5.5  Wave  height,  LF  wave  height  and  setup,  instationary  model  results  vs  measurements  from  LIP  11D, 
test  2E. 
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Figure  5.6  Orbital  velociity  and  long  wave  rms  velocity  ,  instationary  model  results  vs  LIP  1  ID  test  2E. 
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5.5  Dune  erosion  in  large-scale  flume  test. 

The  purpose  of  this  test  was  to  verify  the  dune  erosion  modelling  in  instationary  mode. 

The  model  was  run  for  0.8  hours  of  hydrodynamic  time  with  a  morphological  factor  of  10, 
effectively  representing  a  morphological  simulation  time  of  8  hours. 

A  key  element  in  the  modelling  is  the  avalanching  algorithm;  although  the  surfbeat  waves 
that  are  explicitly  modelled  run  up  and  down  the  upper  beach,  without  a  mechanism  to 
transport  sand  from  dry  to  wet  the  dune  erosion  process  will  not  happen.  A  relatively  simple 
approach,  whereby  an  underwater  critical  slope  of  0.3  and  a  critical  slope  above  water  of  1.0 
were  applied,  proves  to  be  quite  successful  in  representing  the  retreat  of  the  upper  beach  and 
dune  face.  A  grid  resolution  of  1  m  was  applied.  In  Figure  5.7  The  measured  and  modelled 
bed  evolution  is  shown,  which  looks  quite  promising  in  the  upper  region.  The  behaviour  of 
the  bar  at  approx.  135  m  is  not  represented  well;  for  this,  additional  processes  such  as  the 
effect  of  surface  rollers  and  wave  asymmetry /skewness  have  to  be  taken  into  account. 


5.5 

5 

4.5 

4 

3.5 

3 

2.5 


Figure  5.7  Measured  and  modelled  bed  level  after  0,  1,  2,  4  and  8  hours  of  wave  action,  for  a  water  level  of  5.  1 
m  above  the  flume  bottom. 
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6  Conclusions  and  next  steps 


Significant  progress  has  been  made  over  the  last  three  months  in  modelling  the  dune  erosion 
process,  and  a  first  series  of  validation  tests  has  been  carried  out.  The  approach  tested  so  far 
has  been  to  run  in  instationary  mode,  where  wave-group  generated  long  waves  are  generated 
which  represent  the  dominant  motion  in  the  swash  zone.  As  these  motions  are  also  likely  to 
dominate  overwashing  processes  we  are  eager  to  proceed  to  test  cases  where  overwashing 
occurs. 

A  drawback  of  this  approach  may  be  the  required  resolution,  in  the  order  of  metres;  in  ID 
simulations  this  is  no  problem,  in  2DH  simulations  covering  larger  domains  it  may  become 
restrictive,  though  parallellization  can  solve  much  of  this  problem.  Still,  it  is  worth 
considering  alternative  ‘quick-and-dirty’  approaches  based  on  stationary  wave  and  flow 
modelling  combined  with  an  extrapolation  method  for  the  actual  dune  erosion. 

Our  plans  for  the  coming  months  are  summarized  below: 

•  Investigate  SteetzeTs  formulations,  Van  Rijn’s  latest  method 

•  Investigate  original  D3D  “extrapolation  method” 

•  Include  roller  model  and  wave  asymmetry  effects 

•  2D  offshore  absorbing  boundary  conditions  (Ap  and  Jamie) 

•  Include  Q3D  description  of  flow  cf  Reniers  et  al,  2004b 

•  Create  a  collaborative  workspace  on  the  IHE  server. 

•  Start  version  management  through  the  same  server. 

•  Incorporate  all  cases  in  Delft  Hydraulics  test  bed  format 

•  Review  formulations  of  wave  impact  contributions  by  short  waves  and  long  waves 
(e.g.  Overton  and  Fisher). 

•  Future  tests  for  testbed: 

o  New  dune  erosion  test  from  Deltaflume  (available) 
o  Berm  test  (available) 

o  Scheveningen  berm  test  in  Deltaflume  (this  Fall) 
o  Oregon  test  (after  summer) 

o  2D  tests  are  only  in  situ  but  with  little  data.  Possible:  2D  Dauphin, 
Monterey,  East  coast  Florida.  Brad  will  check  with  Abby  Sallenger,  Ap  will 
contact  Brad.  Ad  to  check  with  Thornton,  DJ  to  check  with  CPI.  NCEX  data 
is  also  possibility,  but  little  coastal  damage.  Nobu  Kobayashi  at  Deleware 
University  is  proposing  to  NSF-Intemational  to  conduct  a  lab  study  at  PARI 
Japan  or  at  ORST  Tsunami  tank.  Ap  will  keep  in  contact. 

•  Dano  Roelvink  will  go  to  MORPHOS  meeting  in  Vicksburg  at  November  1-2. 
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